library(dplyr)
library(stringr)
library(multiwayvcov)
library(lmtest)
library(psych)
library(ggplot2)
library(stargazer)
library(margins)
library(vtable)

setwd("E:/Database/Conjoint/Data")
load(file='P4P_Clean.rda')
str(P4P_Clean)

P4P_Clean <- within(P4P_Clean, `Total pay` <- relevel(`Total pay`, ref = "Slightly below average"))
P4P_Clean <- within(P4P_Clean, `Performance bonuses` <- relevel(`Performance bonuses`, ref = "No performance bonuses; fixed salary"))
P4P_Clean <- within(P4P_Clean, `Performance bonuses binary` <- relevel(`Performance bonuses binary`, ref = "No (fixed salary)"))
P4P_Clean <- within(P4P_Clean, `Job performance evaluation` <- relevel(`Job performance evaluation`, ref = "A supervisor evaluation of your work"))
P4P_Clean <- within(P4P_Clean, `Goal based evaluation` <- relevel(`Goal based evaluation`, ref = "No (supervisor based)"))
P4P_Clean <- within(P4P_Clean, `Current community involvement` <- relevel(`Current community involvement`, ref = "Rare participation"))
P4P_Clean <- within(P4P_Clean, `Community income` <- relevel(`Community income`, ref = "Low income"))
P4P_Clean <- within(P4P_Clean, `Community demographics` <- relevel(`Community demographics`, ref = "Mostly white"))
P4P_Clean <- within(P4P_Clean, `Overtime work` <- relevel(`Overtime work`, ref = "Never required"))
P4P_Clean <- within(P4P_Clean, `Key job task` <- relevel(`Key job task`, ref = "Analysis identifying community needs"))
P4P_Clean <- within(P4P_Clean, Race_B <- relevel(Race_B, ref = "White"))
P4P_Clean <- within(P4P_Clean, Race_C <- relevel(Race_C, ref = "White"))
P4P_Clean <- within(P4P_Clean, Gender <- relevel(Gender, ref = "Male"))


# 1. Descriptive Statistics (Table 2)
P4P_Clean_Des <- P4P_Clean[c(1:1501),] %>%
  mutate(Education=as.factor(Education),
         Employment=as.factor(Employment),
         Income=as.factor(Income))

vtable::st(P4P_Clean_Des, 
           vars = c('Gender','Race_C','Education','Employment','Income','Party'), 
           digits = 3)
vtable::st(P4P_Clean_Des, 
           vars = c('Gender','Race_C','Education','Employment','Income','Party'),
           group='Study.ID', digits = 3)

psych::describe(P4P_Clean_Des$Age)
psych::describeBy(P4P_Clean_Des$Age,P4P_Clean_Des$Study.ID)


# 2. Regression Analysis
# 2.1. All Attributes
PFP1 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task`, 
           data=P4P_Clean)
summary(PFP1)
PFP1.cl <- cluster.vcov(PFP1, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(PFP1, PFP1.cl)
Twoway.se <- sqrt(diag(PFP1.cl))


# 2.1.1 Levels of P4P
PFP2 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task`, 
           data=P4P_Clean)
summary(PFP2)
PFP2.cl <- cluster.vcov(PFP2, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(PFP2, PFP2.cl)
Twoway.se.NotCollapse <- sqrt(diag(PFP2.cl))


# 2.1.2 Interaction w/ Goal based evaluation
PFP3 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + `Performance bonuses binary`:`Goal based evaluation`, 
           data=P4P_Clean)
summary(PFP3)
PFP3.cl <- cluster.vcov(PFP3, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(PFP3, PFP3.cl)
Twoway.se.PG <- sqrt(diag(PFP3.cl))


# 2.1.3 Interaction w/ Risk Attitude
PFP4 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + GeneralRiskB + `Performance bonuses binary`:GeneralRiskB, 
           data=P4P_Clean)
summary(PFP4)
PFP4.cl <- cluster.vcov(PFP4, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(PFP4, PFP4.cl)
Twoway.se.PRG <- sqrt(diag(PFP4.cl))

PFP5 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + LotteryRiskB + `Performance bonuses binary`:LotteryRiskB, 
           data=P4P_Clean)
summary(PFP5)
PFP5.cl <- cluster.vcov(PFP5, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(PFP5, PFP5.cl)
Twoway.se.PRL <- sqrt(diag(PFP5.cl))

## Report (Table A1)
stargazer(PFP1, PFP2,PFP3,PFP4,PFP5,
          se=list(Twoway.se,Twoway.se.NotCollapse,Twoway.se.PG,Twoway.se.PRG,Twoway.se.PRL), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%)",
                             "Performance bonuses: moderate (10%)",
                             "Performance bonuses: small (5%)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "(A) X (B)",
                             "General Risk: High Risk", 
                             "(A) X General Risk: High Risk",
                             "Lottery Risk: High Risk",
                             "(A) X Lottery Risk: High Risk"),
          out="Main_P4P.txt") 


# 2.2 Race
table(P4P_Clean$Race_B,exclude=NULL)
table(P4P_Clean$Race_C,exclude=NULL)

P4P_Data_Race <- P4P_Clean %>%
  subset(Race_B!="Prefer not to answer")

Race1 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
              `Current community involvement` + `Community income` + `Community demographics` +
              `Overtime work` + `Key job task` + Race_B + `Performance bonuses binary`*Race_B, 
            data=P4P_Data_Race)
summary(Race1)
Race1.cl <- cluster.vcov(Race1, cbind(P4P_Data_Race$ID,P4P_Data_Race$Treatment.ID))
coeftest(Race1, Race1.cl)
Twoway.se.RaceB <- sqrt(diag(Race1.cl))


# 2.2.1 Levels of P4P
Race2 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
              `Current community involvement` + `Community income` + `Community demographics` +
              `Overtime work` + `Key job task` + Race_B + `Performance bonuses`:Race_B, 
            data=P4P_Data_Race)
summary(Race2)
Race2.cl <- cluster.vcov(Race2, cbind(P4P_Data_Race$ID, P4P_Data_Race$Treatment.ID))
coeftest(Race2, Race2.cl)
Twoway.se.RaceB.Prop <- sqrt(diag(Race2.cl))

# 2.2.2 Interaction w/ Goal based evaluation
Race3 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
              `Current community involvement` + `Community income` + `Community demographics` +
              `Overtime work` + `Key job task` + Race_B + 
              `Performance bonuses binary`:Race_B + `Goal based evaluation`:Race_B +
              `Performance bonuses binary`:`Goal based evaluation` +
              `Performance bonuses binary`:`Goal based evaluation`:Race_B, 
            data=P4P_Data_Race)
summary(Race3)
Race3.cl <- cluster.vcov(Race3, cbind(P4P_Data_Race$ID, P4P_Data_Race$Treatment.ID))
coeftest(Race3, Race3.cl)
Twoway.se.EVAL <- sqrt(diag(Race3.cl))

# 2.2.3 Interaction w/ Risk Attitude
Race4 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
              `Current community involvement` + `Community income` + `Community demographics` +
              `Overtime work` + `Key job task` + Race_B + GeneralRiskB +
              `Performance bonuses binary`*Race_B + `Performance bonuses binary`*GeneralRiskB, 
            data=P4P_Data_Race)
summary(Race4)
Race4.cl <- cluster.vcov(Race4, cbind(P4P_Data_Race$ID, P4P_Data_Race$Treatment.ID))
coeftest(Race4, Race4.cl)
Twoway.se.RaceB.RiskG <- sqrt(diag(Race4.cl))

Race5 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
              `Current community involvement` + `Community income` + `Community demographics` +
              `Overtime work` + `Key job task` + Race_B + LotteryRiskB +
              `Performance bonuses binary`*Race_B + `Performance bonuses binary`*LotteryRiskB, 
            data=P4P_Data_Race)
summary(Race5)
Race5.cl <- cluster.vcov(Race5, cbind(P4P_Data_Race$ID, P4P_Data_Race$Treatment.ID))
coeftest(Race5, Race5.cl)
Twoway.se.RaceB.RiskL <- sqrt(diag(Race5.cl))


## Report (Table A2)
stargazer(Race1,Race2,Race3,Race4,Race5,
          se=list(Twoway.se.RaceB,Twoway.se.RaceB.Prop,Twoway.se.EVAL,Twoway.se.RaceB.RiskG,Twoway.se.RaceB.RiskL), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%) (A1)",
                             "Performance bonuses: moderate (10%) (A2)",
                             "Performance bonuses: small (5%) (A3)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "Race - binary: Minority (C)",
                             "General Risk: High Risk",
                             "Lottery Risk: High Risk",
                             "(A) X (C)",
                             "(A1) X (C)",
                             "(A2) X (C)",
                             "(A3) X (C)",
                             "(B) X (C)",
                             "(A) X (B)",
                             "(A) X (B) X (C)",
                             "(A) X General Risk: High Risk",
                             "(A) X Lottery Risk: High Risk"),
          out="Race_Main_P4P.txt") 


# 2.2.4 Race - Specific
Race6 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
              `Current community involvement` + `Community income` + `Community demographics` +
              `Overtime work` + `Key job task` + Race_C + `Performance bonuses binary`:Race_C, 
            data=P4P_Data_Race)
summary(Race6)
Race6.cl <- cluster.vcov(Race6, cbind(P4P_Data_Race$ID, P4P_Data_Race$Treatment.ID))
coeftest(Race6, Race6.cl)
Twoway.se.RaceC <- sqrt(diag(Race6.cl))

Race7 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
              `Current community involvement` + `Community income` + `Community demographics` +
              `Overtime work` + `Key job task` + Race_C + `Performance bonuses`:Race_C, 
            data=P4P_Data_Race)
summary(Race7)
Race7.cl <- cluster.vcov(Race7, cbind(P4P_Data_Race$ID, P4P_Data_Race$Treatment.ID))
coeftest(Race7, Race7.cl)
Twoway.se.RaceC.Prop <- sqrt(diag(Race7.cl))

## Report (Table A3)
stargazer(Race6, Race7,
          se=list(Twoway.se.RaceC,Twoway.se.RaceC.Prop), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%) (A1)",
                             "Performance bonuses: moderate (10%) (A2)",
                             "Performance bonuses: small (5%) (A3)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "Race - Asian", "Race - Black", "Race - Hispanic", "Race - Other",
                             "(A) X Race - Asian", "(A) X Race - Black", "(A) X Race - Hispanic", "(A) X Race - Other",
                             "(A1) X Race - Asian", "(A2) X Race - Asian", "(A3) X Race - Asian",
                             "(A1) X Race - Black", "(A2) X Race - Black", "(A3) X Race - Black",
                             "(A1) X Race - Hispanic", "(A2) X Race - Hispanic", "(A3) X Race - Hispanic",
                             "(A1) X Race - Other", "(A2) X Race - Other", "(A3) X Race - Other"),
          out="Race_Separate_P4P.txt") 


# 2.3 Gender
P4P_Data_Gender <- P4P_Clean %>%
  subset(Gender=="Female" | Gender=="Male") %>%
  mutate(Gender=factor(Gender))

Gen1 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + `Performance bonuses binary`:Gender, 
           data=P4P_Data_Gender)
summary(Gen1)
Gen1.cl <- cluster.vcov(Gen1, cbind(P4P_Data_Gender$ID, P4P_Data_Gender$Treatment.ID))
coeftest(Gen1, Gen1.cl)
Twoway.se.Gender <- sqrt(diag(Gen1.cl))

# 2.3.1 Levels of P4P
Gen2 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + `Performance bonuses`:Gender, 
           data=P4P_Data_Gender)
summary(Gen2)
Gen2.cl <- cluster.vcov(Gen2, cbind(P4P_Data_Gender$ID, P4P_Data_Gender$Treatment.ID))
coeftest(Gen2, Gen2.cl)
Twoway.se.Gender.Prop <- sqrt(diag(Gen2.cl))

# 2.3.2 Interaction w/ Goal based evaluation
Gen3 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + `Performance bonuses binary`*`Goal based evaluation`*Gender, 
           data=P4P_Data_Gender)
summary(Gen3)
Gen3.cl <- cluster.vcov(Gen3, cbind(P4P_Data_Gender$ID, P4P_Data_Gender$Treatment.ID))
coeftest(Gen3, Gen3.cl)
Twoway.se.EVAL.Gender <- sqrt(diag(Gen3.cl))

# 2.3.3 Interaction w/ Risk Attitude
Gen4 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + GeneralRiskB +
             `Performance bonuses binary`:Gender + `Performance bonuses binary`:GeneralRiskB, 
           data=P4P_Data_Gender)
summary(Gen4)
Gen4.cl <- cluster.vcov(Gen4, cbind(P4P_Data_Gender$ID, P4P_Data_Gender$Treatment.ID))
coeftest(Gen4, Gen4.cl)
Twoway.se.Gender.RiskG <- sqrt(diag(Gen4.cl))

Gen5 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + LotteryRiskB +
             `Performance bonuses binary`:Gender + `Performance bonuses binary`:LotteryRiskB, 
           data=P4P_Data_Gender)
summary(Gen5)
Gen5.cl <- cluster.vcov(Gen5, cbind(P4P_Data_Gender$ID, P4P_Data_Gender$Treatment.ID))
coeftest(Gen5, Gen5.cl)
Twoway.se.Gender.RiskL <- sqrt(diag(Gen5.cl))


## Report (Table A4)
stargazer(Gen1,Gen2,Gen3,Gen4,Gen5,
          se=list(Twoway.se.Gender,Twoway.se.Gender.Prop,Twoway.se.EVAL.Gender,Twoway.se.Gender.RiskG,Twoway.se.Gender.RiskL), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%) (A1)",
                             "Performance bonuses: moderate (10%) (A2)",
                             "Performance bonuses: small (5%) (A3)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "Gender: Women (C)",
                             "(A) X (B)",
                             "General Risk: High Risk", "Lottery Risk: High Risk",
                             "(A) X (C)",
                             "(A1) X (C)",
                             "(A2) X (C)",
                             "(A3) X (C)",
                             "(B) X (C)",
                             "(A) X (B) X (C)",
                             "(A) X General Risk: High Risk", "(A) X Lottery Risk: High Risk"),
          out="Gender_Main_P4P.txt") 


# 2.4 Age
Age1 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + `Performance bonuses binary`:Age_B, 
           data=P4P_Clean)
summary(Age1)
Age1.cl <- cluster.vcov(Age1, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(Age1, Age1.cl)
Twoway.se.Age <- sqrt(diag(Age1.cl))

# 2.4.1 Levels of P4P
Age2 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + `Performance bonuses`:Age_B, 
           data=P4P_Clean)
summary(Age2)
Age2.cl <- cluster.vcov(Age2, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(Age2, Age2.cl)
Twoway.se.Age.Prop <- sqrt(diag(Age2.cl))

# 2.4.2 Interaction w/ Goal based evaluation
Age3 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + `Performance bonuses binary`*`Goal based evaluation`*Age_B, 
           data=P4P_Clean)
summary(Age3)
Age3.cl <- cluster.vcov(Age3, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(Age3, Age3.cl)
Twoway.se.EVAL.Age <- sqrt(diag(Age3.cl))

# 2.4.3 Interaction w/ Risk Attitude
Age4 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + GeneralRiskB +
             `Performance bonuses binary`:Age_B + `Performance bonuses binary`:GeneralRiskB, 
           data=P4P_Clean)
summary(Age4)
Age4.cl <- cluster.vcov(Age4, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(Age4, Age4.cl)
Twoway.se.Age.RiskG <- sqrt(diag(Age4.cl))

Age5 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + LotteryRiskB +
             `Performance bonuses binary`:Age_B + `Performance bonuses binary`:LotteryRiskB, 
           data=P4P_Clean)
summary(Age5)
Age5.cl <- cluster.vcov(Age5, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(Age5, Age5.cl)
Twoway.se.Age.RiskL <- sqrt(diag(Age5.cl))

## Report (Table A5)
stargazer(Age1,Age2,Age3,Age4,Age5,
          se=list(Twoway.se.Age,Twoway.se.Age.Prop,Twoway.se.EVAL.Age,Twoway.se.Age.RiskG,Twoway.se.Age.RiskL), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%) (A1)",
                             "Performance bonuses: moderate (10%) (A2)",
                             "Performance bonuses: small (5%) (A3)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "Age: Older (C)",
                             "(A) X (B)",
                             "General Risk: High Risk", "Lottery Risk: High Risk",
                             "(A) X (C)",
                             "(A1) X (C)",
                             "(A2) X (C)",
                             "(A3) X (C)",
                             "(B) X (C)",
                             "(A) X (B) X (C)",
                             "(A) X General Risk: High Risk", "(A) X Lottery Risk: High Risk"),
          out="Age_Main_P4P.txt") 


# 3. Plot
# 3.1 Overall
Overall.margin <- margins::margins(PFP1, vcov = cluster.vcov(PFP1, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))) %>%
  summary() %>%
  mutate(Value=c("Mostly African American","Mostly Hispanic","Multiracial",
                 "Average income","High income",
                 "Frequent participation","Moderate participation",
                 "Yes (Goal based)",
                 "Coordination with community groups and organizations",
                 "Direct interaction with community residents","Teamwork with peers and supervisors",
                 "Frequently required","Occasionally required",
                 "Yes (bonuses)",
                 "About average","Slightly above average")) %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Base.Factor <- data.frame(
  Value=c("Total pay:", "(Baseline=Slightly below average)","About average","Slightly above average",
          "Performance bonuses binary:", "(Baseline=No (Fixed salary))","Yes (bonuses)",
          "Goal based evaluation:","(Baseline=No (Supervisor based)","Yes (Goal based)",
          "Current community involvement:","(Baseline=Rare participation","Frequent participation","Moderate participation",
          "Community income:","(Baseline=Low income)","Average income","High income",
          "Community demographics:","(Baseline=Mostly white)","Mostly African American","Mostly Hispanic","Multiracial",
          "Overtime work:","(Baseline=Never required)","Frequently required","Occasionally required",
          "Key job task:","(Baseline=Analysis identifying community needs)","Coordination with community groups and organizations",
          "Direct interaction with community residents","Teamwork with peers and supervisors"))

Overall.table <- merge(x=Base.Factor, y=Overall.margin,
                       by='Value', all.x=TRUE)
Overall.table$Value <- factor(Overall.table$Value,
                              levels=c("Total pay:", "(Baseline=Slightly below average)","About average","Slightly above average",
                                       "Performance bonuses binary:", "(Baseline=No (Fixed salary))","Yes (bonuses)",
                                       "Goal based evaluation:","(Baseline=No (Supervisor based)","Yes (Goal based)",
                                       "Current community involvement:","(Baseline=Rare participation","Frequent participation","Moderate participation",
                                       "Community income:","(Baseline=Low income)","Average income","High income",
                                       "Community demographics:","(Baseline=Mostly white)","Mostly African American","Mostly Hispanic","Multiracial",
                                       "Overtime work:","(Baseline=Never required)","Frequently required","Occasionally required",
                                       "Key job task:","(Baseline=Analysis identifying community needs)","Coordination with community groups and organizations",
                                       "Direct interaction with community residents","Teamwork with peers and supervisors"))

# Figure 2
ggplot(Overall.table, aes(x = AME, y = Value)) + 
  geom_vline(aes(xintercept = 0), size = .25, linetype = "dashed") + 
  geom_errorbarh(data = Overall.table, aes(xmax = upper, xmin = lower), size = .3, height = .6, position = position_nudge(y = 0), color="black") +
  geom_errorbarh(data = Overall.table, aes(xmax = upper_90, xmin = lower_90), size = .99, height = .6, position = position_nudge(y = 0), color="black") +
  geom_point(data = Overall.table, shape=19, size = 3, color='black') +
  scale_y_discrete(limits=rev) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(axis.text.y = element_text(face = 
                                     c("plain","plain","plain","plain","bold","plain","plain","plain","bold","plain","plain","plain","plain","bold",
                                       "plain","plain","plain","bold","plain","plain","plain","bold","plain","plain","bold","plain","plain","bold",
                                       "plain","plain","plain","bold"))) +
  ylab("") +
  xlab("AMCE: Effect on Pr(Job Preferred)") +
  theme(axis.text.y = element_text(size=12, colour = "black")) +
  theme(text = element_text(size=12, colour = "black"))


# 3.2 Interactions
Interaction.Race <- margins::margins(Race1, 
                                     vcov = cluster.vcov(Race1, cbind(P4P_Data_Race$ID, P4P_Data_Race$Treatment.ID)), 
                                     variables = "Performance bonuses binary",
                                     at=list(Race_B=c("White","Minority"))) %>%
  summary() %>% rename(Value=Race_B) %>% mutate(group="Race (binary)") %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Interaction.RaceC <- margins::margins(Race6, 
                                      vcov = cluster.vcov(Race6, cbind(P4P_Data_Race$ID, P4P_Data_Race$Treatment.ID)), 
                                      variables = "Performance bonuses binary",
                                      at=list(Race_C=c("White","Black","Hispanic","Asian","Other"))) %>%
  summary() %>% rename(Value=Race_C) %>% mutate(group="Race") %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Interaction.Gender <- margins::margins(Gen1, 
                                       vcov = cluster.vcov(Gen1, cbind(P4P_Data_Gender$ID, P4P_Data_Gender$Treatment.ID)), 
                                       variables = "Performance bonuses binary",
                                       at=list(Gender=c("Female","Male"))) %>%
  summary() %>% rename(Value=Gender) %>% mutate(group="Gender") %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Interaction.Age <- margins::margins(Age1, 
                                    vcov = cluster.vcov(Age1, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID)), 
                                    variables = "Performance bonuses binary",
                                    at=list(Age_B=c(1,0))) %>%
  summary() %>% mutate(Value=ifelse(Age_B==1,"Older","Younger")) %>% 
  mutate(group="Age") %>%
  dplyr::select(factor,Value,AME,SE,z,p,lower,upper,group) %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Interaction <- rbind(Interaction.Race,Interaction.RaceC,Interaction.Gender,Interaction.Age)

# Figure 3
ggplot(Interaction, aes(x = AME, y = group)) + 
  geom_vline(aes(xintercept = 0), size = .25, linetype = "dashed") + 
  geom_errorbarh(data = filter(Interaction, Value=="White" & group=="Race (binary)"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="White" & group=="Race (binary)"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_point(data = filter(Interaction, Value== "White" & group=="Race (binary)"), shape=19, size = 3, position = position_nudge(y = 0.1), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="Minority"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="Minority"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_point(data = filter(Interaction, Value== "Minority"), shape=19, size = 3, position = position_nudge(y = -0.1), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="White" & group=="Race"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.4), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="White" & group=="Race"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.4), color="black") +
  geom_point(data = filter(Interaction, Value== "White" & group=="Race"), shape=19, size = 3, position = position_nudge(y = 0.4), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="Asian"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.2), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="Asian"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.2), color="black") +
  geom_point(data = filter(Interaction, Value== "Asian"), shape=19, size = 3, position = position_nudge(y = 0.2), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="Black"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.0), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="Black"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.0), color="black") +
  geom_point(data = filter(Interaction, Value== "Black"), shape=19, size = 3, position = position_nudge(y = 0.0), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="Hispanic"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.2), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="Hispanic"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.2), color="black") +
  geom_point(data = filter(Interaction, Value== "Hispanic"), shape=19, size = 3, position = position_nudge(y = -0.2), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="Other"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.4), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="Other"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.4), color="black") +
  geom_point(data = filter(Interaction, Value== "Other"), shape=19, size = 3, position = position_nudge(y = -0.4), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="Male"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="Male"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_point(data = filter(Interaction, Value== "Male"), shape=19, size = 3, position = position_nudge(y = 0.1), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="Female"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="Female"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_point(data = filter(Interaction, Value== "Female"), shape=19, size = 3, position = position_nudge(y = -0.1), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="Older"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="Older"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_point(data = filter(Interaction, Value== "Older"), shape=19, size = 3, position = position_nudge(y = 0.1), color='black') +
  geom_errorbarh(data = filter(Interaction, Value=="Younger"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_errorbarh(data = filter(Interaction, Value=="Younger"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_point(data = filter(Interaction, Value== "Younger"), shape=19, size = 3, position = position_nudge(y = -0.1), color='black') +
  coord_cartesian(xlim = c(-0.2, 0.2)) +
  ylab("") +
  xlab("AMCE: Effect of Pay-for-Performance on Pr(Job Preferred)") +
  annotate(geom="text", x=-0.18, y=4.1, label="White", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=3.9, label="Minority", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=3.4, label="White", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=3.2, label="Asian", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=3.0, label="African American", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=2.8, label="Hispanic", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=2.6, label="Other", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=2.1, label="Men", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=1.9, label="Women", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=1.1, label="Older", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.18, y=0.9, label="Younger", colour="black", fontface = 'italic', size=4) +
  theme(axis.text.y=element_text(size=12,face="bold",colour="black"))


## Other Interactions (Appendix Figure A1)
Evaluation <- margins::margins(PFP3, 
                               vcov = cluster.vcov(PFP3, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID)), 
                               variables = "Performance bonuses binary",
                               at=list(`Goal based evaluation`=c("Yes (goal based)","No (supervisor based)"))) %>%
  summary() %>% mutate(Value=ifelse(`Goal based evaluation`=="Yes (goal based)","Goal based evaluation","Supervisor based evaluation")) %>% 
  mutate(group="Performance Evaluation") %>%
  dplyr::select(factor,Value,AME,SE,z,p,lower,upper,group) %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Risk1 <- margins::margins(PFP4, 
                          vcov = cluster.vcov(PFP4, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID)), 
                          variables = "Performance bonuses binary",
                          at=list(GeneralRiskB=c(1,0))) %>%
  summary() %>% mutate(Value=ifelse(GeneralRiskB==1,"High","Low")) %>% 
  mutate(group="General Risk") %>%
  dplyr::select(factor,Value,AME,SE,z,p,lower,upper,group) %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Risk2 <- margins::margins(PFP5, 
                          vcov = cluster.vcov(PFP5, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID)), 
                          variables = "Performance bonuses binary",
                          at=list(LotteryRiskB=c(1,0))) %>%
  summary() %>% mutate(Value=ifelse(LotteryRiskB==1,"High","Low")) %>% 
  mutate(group="Lottery Risk") %>%
  dplyr::select(factor,Value,AME,SE,z,p,lower,upper,group) %>%
  mutate(upper_90=AME+qnorm(0.95)*SE, lower_90=AME-qnorm(0.95)*SE)

Interaction_Other <- rbind(Evaluation,Risk1,Risk2)
level_order <- c("Lottery Risk","General Risk","Performance Evaluation")

# Figure A2
ggplot(Interaction_Other, aes(x = AME, y = group)) + 
  geom_vline(aes(xintercept = 0), size = .25, linetype = "dashed") + 
  geom_errorbarh(data = filter(Interaction_Other, Value=="Goal based evaluation"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_errorbarh(data = filter(Interaction_Other, Value=="Goal based evaluation"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_point(data = filter(Interaction_Other, Value== "Goal based evaluation"), shape=19, size = 3, position = position_nudge(y = 0.1), color='black') +
  geom_errorbarh(data = filter(Interaction_Other, Value=="Supervisor based evaluation"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_errorbarh(data = filter(Interaction_Other, Value=="Supervisor based evaluation"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_point(data = filter(Interaction_Other, Value== "Supervisor based evaluation"), shape=19, size = 3, position = position_nudge(y = -0.1), color='black') +
  geom_errorbarh(data = filter(Interaction_Other, Value=="High" & group=="General Risk"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_errorbarh(data = filter(Interaction_Other, Value=="High" & group=="General Risk"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_point(data = filter(Interaction_Other, Value== "High" & group=="General Risk"), shape=19, size = 3, position = position_nudge(y = 0.1), color='black') +
  geom_errorbarh(data = filter(Interaction_Other, Value=="Low" & group=="General Risk"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_errorbarh(data = filter(Interaction_Other, Value=="Low" & group=="General Risk"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_point(data = filter(Interaction_Other, Value== "Low" & group=="General Risk"), shape=19, size = 3, position = position_nudge(y = -0.1), color='black') +
  geom_errorbarh(data = filter(Interaction_Other, Value=="High" & group=="Lottery Risk"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_errorbarh(data = filter(Interaction_Other, Value=="High" & group=="Lottery Risk"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.1), color="black") +
  geom_point(data = filter(Interaction_Other, Value== "High" & group=="Lottery Risk"), shape=19, size = 3, position = position_nudge(y = 0.1), color='black') +
  geom_errorbarh(data = filter(Interaction_Other, Value=="Low" & group=="Lottery Risk"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_errorbarh(data = filter(Interaction_Other, Value=="Low" & group=="Lottery Risk"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.1), color="black") +
  geom_point(data = filter(Interaction_Other, Value== "Low" & group=="Lottery Risk"), shape=19, size = 3, position = position_nudge(y = -0.1), color='black') +
  coord_cartesian(xlim = c(-0.15, 0.15)) +
  scale_y_discrete(limits = level_order) +
  ylab("") +
  xlab("AMCE: Effect of Pay-for-Performance on Pr(Job Preferred)") +
  annotate(geom="text", x=-0.11, y=3.2, label="Goal based evaluation", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.11, y=2.8, label="Supervisor based evaluation", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.14, y=2.1, label="High", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.14, y=1.9, label="Low", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.14, y=1.1, label="High", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.14, y=0.9, label="Low", colour="black", fontface = 'italic', size=4) +
  theme(axis.text.y=element_text(size=12,face="bold",colour="black"))


## w/ Reference (Appendix Figure A2)
RaceB_coef <- data.frame(coeftest(Race1, Race1.cl)[19,1])
names(RaceB_coef) <- "Coef"
RaceB_se <- data.frame(coeftest(Race1, Race1.cl)[19,2])
names(RaceB_se) <- "SE"
RaceB_coef_se <- cbind(RaceB_coef,RaceB_se)

RaceC_coef <- data.frame(coeftest(Race6, Race6.cl)[22:25,1])
names(RaceC_coef) <- "Coef"
RaceC_se <- data.frame(coeftest(Race6, Race6.cl)[22:25,2])
names(RaceC_se) <- "SE"
RaceC_coef_se <- cbind(RaceC_coef,RaceC_se)

Gender_coef <- data.frame(coeftest(Gen1, Gen1.cl)[19,1])
names(Gender_coef) <- "Coef"
Gender_se <- data.frame(coeftest(Gen1, Gen1.cl)[19,2])
names(Gender_se) <- "SE"
Gender_coef_se <- cbind(Gender_coef,Gender_se)

Age_coef <- data.frame(coeftest(Age1, Age1.cl)[19,1])
names(Age_coef) <- "Coef"
Age_se <- data.frame(coeftest(Age1, Age1.cl)[19,2])
names(Age_se) <- "SE"
Age_coef_se <- cbind(Age_coef,Age_se)

Int_Ref <-rbind(RaceB_coef_se,RaceC_coef_se,Gender_coef_se,Age_coef_se) %>%
  mutate(upper=Coef+qnorm(0.975)*SE, lower=Coef-qnorm(0.975)*SE) %>%
  mutate(upper_90=Coef+qnorm(0.95)*SE, lower_90=Coef-qnorm(0.95)*SE) %>%
  mutate(Value=c("Minority (ref. White)",
                 "Asian (ref. White)","African American (ref. White)","Hispanic (ref. White)","Other (ref. White)",
                 "Women (ref. Men)",
                 "Older (ref. Younger)")) %>%
  mutate(group=c("Race (binary)",
                 "Race","Race","Race","Race",
                 "Gender",
                 "Age"))

# Figure A2
ggplot(Int_Ref, aes(x = Coef, y = group)) + 
  geom_vline(aes(xintercept = 0), size = .25, linetype = "dashed") + 
  geom_errorbarh(data = filter(Int_Ref, Value=="Minority (ref. White)" & group=="Race (binary)"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.0), color="black") +
  geom_errorbarh(data = filter(Int_Ref, Value=="Minority (ref. White)" & group=="Race (binary)"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.0), color="black") +
  geom_point(data = filter(Int_Ref, Value== "Minority (ref. White)" & group=="Race (binary)"), shape=19, size = 3, position = position_nudge(y = 0.0), color='black') +
  geom_errorbarh(data = filter(Int_Ref, Value=="Asian (ref. White)"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.45), color="black") +
  geom_errorbarh(data = filter(Int_Ref, Value=="Asian (ref. White)"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.45), color="black") +
  geom_point(data = filter(Int_Ref, Value== "Asian (ref. White)"), shape=19, size = 3, position = position_nudge(y = 0.45), color='black') +
  geom_errorbarh(data = filter(Int_Ref, Value=="African American (ref. White)"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = 0.15), color="black") +
  geom_errorbarh(data = filter(Int_Ref, Value=="African American (ref. White)"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = 0.15), color="black") +
  geom_point(data = filter(Int_Ref, Value== "African American (ref. White)"), shape=19, size = 3, position = position_nudge(y = 0.15), color='black') +
  geom_errorbarh(data = filter(Int_Ref, Value=="Hispanic (ref. White)"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.15), color="black") +
  geom_errorbarh(data = filter(Int_Ref, Value=="Hispanic (ref. White)"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.15), color="black") +
  geom_point(data = filter(Int_Ref, Value== "Hispanic (ref. White)"), shape=19, size = 3, position = position_nudge(y = -0.15), color='black') +
  geom_errorbarh(data = filter(Int_Ref, Value=="Other (ref. White)"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.45), color="black") +
  geom_errorbarh(data = filter(Int_Ref, Value=="Other (ref. White)"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.45), color="black") +
  geom_point(data = filter(Int_Ref, Value== "Other (ref. White)"), shape=19, size = 3, position = position_nudge(y = -0.45), color='black') +
  geom_errorbarh(data = filter(Int_Ref, Value=="Women (ref. Men)"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.0), color="black") +
  geom_errorbarh(data = filter(Int_Ref, Value=="Women (ref. Men)"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.0), color="black") +
  geom_point(data = filter(Int_Ref, Value== "Women (ref. Men)"), shape=19, size = 3, position = position_nudge(y = -0.0), color='black') +
  geom_errorbarh(data = filter(Int_Ref, Value=="Older (ref. Younger)"), aes(xmax = upper, xmin = lower, col=group), size = .3, height = .1, position = position_nudge(y = -0.0), color="black") +
  geom_errorbarh(data = filter(Int_Ref, Value=="Older (ref. Younger)"), aes(xmax = upper_90, xmin = lower_90, col=group), size = .99, height = .1, position = position_nudge(y = -0.0), color="black") +
  geom_point(data = filter(Int_Ref, Value== "Older (ref. Younger)"), shape=19, size = 3, position = position_nudge(y = -0.0), color='black') +
  coord_cartesian(xlim = c(-0.2, 0.2)) +
  ylab("") +
  xlab("AMCE: Difference in Effect of Pay-for-Performance on Pr(Job Preferred)") +
  annotate(geom="text", x=-0.16, y=4.0, label="Minority (vs. White)", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.16, y=3.45, label="Asian (vs. White)", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.16, y=3.15, label="African American (vs. White)", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.16, y=2.85, label="Hispanic (vs. White)", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.16, y=2.55, label="Other (vs. White)", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.16, y=2.0, label="Women (ref. Men)", colour="black", fontface = 'italic', size=4) +
  annotate(geom="text", x=-0.16, y=1.0, label="Older (ref. Younger)", colour="black", fontface = 'italic', size=4) +
  theme(axis.text.y=element_text(size=12,face="bold",colour="black"))


# 4. Assumptions 
# 4.1 Contest No
PFP_INT_Contest <- lm(Chosen_Job ~ `Total pay`*Contest_no + `Performance bonuses binary`*Contest_no + `Goal based evaluation`*Contest_no +
                        `Current community involvement`*Contest_no + `Community income`*Contest_no + `Community demographics`*Contest_no +
                        `Overtime work`*Contest_no + `Key job task`*Contest_no, 
                      data=P4P_Clean)
summary(PFP_INT_Contest)
PFP_INT_Contest.cl <- cluster.vcov(PFP_INT_Contest, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(PFP_INT_Contest, PFP_INT_Contest.cl)
Twoway.se.Contest <- sqrt(diag(PFP_INT_Contest.cl))

waldtest(PFP_INT_Contest, PFP1, vcov = PFP_INT_Contest.cl)
#### NOTE ####
# As the above test does not consider the number of clusters,
# we calculated the significance of F-statistics using Stata.
# F(34,126) = 1.50, Prob > F = 0.0571
# See Baum, C. F., Nichols, A., & Schaffer, M. E. (2010, September). 
# Evaluating one-way and two-way cluster-robust covariance matrix estimates. 
# In BOS10 Stata Conference (Vol. 11). Stata Users Group.
# "When a cluster-robust VCE has been calculated, 
# Wald t or F test statistics should take account of the number of clusters,
# rather than relying on the asymptotically behavior of the statistic as N → ∞. 
# The approach that Stata follows involves using the t distribution with 
# G − 1 degrees of freedom rather than N − k degrees of freedom" (p.22)

# 4.2 Profile No
PFP_INT_Profile <- lm(Chosen_Job ~ `Total pay`*Profile_no + `Performance bonuses binary`*Profile_no + `Goal based evaluation`*Profile_no +
                        `Current community involvement`*Profile_no + `Community income`*Profile_no + `Community demographics`*Profile_no +
                        `Overtime work`*Profile_no + `Key job task`*Profile_no, 
                      data=P4P_Clean)
summary(PFP_INT_Profile)
PFP_INT_Profile.cl <- cluster.vcov(PFP_INT_Profile, cbind(P4P_Clean$ID, P4P_Clean$Treatment.ID))
coeftest(PFP_INT_Profile, PFP_INT_Profile.cl)
Twoway.se.Profile <- sqrt(diag(PFP_INT_Profile.cl))

waldtest(PFP_INT_Profile, PFP1, vcov = PFP_INT_Profile.cl)
# See the NOTE above #
# F(16,126) = 1.26, Prob > F = 0.2323

## Report (Appendix Table A6 and A7)
stargazer(PFP_INT_Contest,
          se=list(Twoway.se.Contest), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          out="Contest_P4P.txt") 

stargazer(PFP_INT_Profile,
          se=list(Twoway.se.Profile), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          out="Profile_P4P.txt") 


# 4.3 Randomization (Appendix Figure A3 & A4)
groups <- c(quo(`Total pay`), quo(`Performance bonuses`), quo(`Job performance evaluation`),
            quo(`Current community involvement`),quo(`Community income`),quo(`Community demographics`),
            quo(`Overtime work`),quo(`Key job task`))

Level_Randomization = list()

for (i in seq_along(groups)) {
  N <- P4P_Clean %>% 
    group_by(!!groups[[i]]) %>% 
    summarise(N = n()) %>%
    rename(Value=!!groups[[i]]) %>%
    mutate(Attr=groups[i])
  Level_Randomization[[i]] <- N
}

Level_Randomization_N <- do.call(rbind,Level_Randomization)
Level_Randomization_Prop <- Level_Randomization_N %>% mutate(Prop=N/9006)

Level_Randomization_Label <- data.frame(trimws(gsub("[[:punct:]]","",Level_Randomization_N$Attr),"both"))
names(Level_Randomization_Label) <- "Attribute"

Short.Level <- data.frame(c("Below average","Average","Above average",
                            "Fixed","20%","10%","5%",
                            "Supervisor","Attendance","Changes","Satisfaction",
                            "Rare","Frequent","Moderate",
                            "Low","Average","High",                                
                            "White","Black","Hispanic","Multiracial",
                            "Never","Frequently","Occasionally",
                            "Analysis","Coordination","Interaction","Teamwork"))
names(Short.Level) <- "Short.level"
Level_Randomization_Final <- cbind(Level_Randomization_Prop,Short.Level,Level_Randomization_Label)

ggplot(Level_Randomization_Final, aes(x = Short.level, y = Prop)) +
  geom_col(size=1) +
  facet_wrap(~ Attribute, scales = "free") +
  ylab("Percentage") +
  xlab("Value") +
  theme(axis.text.x = element_text(angle = 8, vjust = 0.5, hjust=0.5))


# Treatment ID
P4P_Data_TID <- P4P_Clean %>% group_by(Treatment.ID) %>%
  slice_max(ID, with_ties = TRUE)

for (i in seq_along(groups)) {
  N <- P4P_Data_TID %>% 
    group_by(!!groups[[i]]) %>% 
    summarise(N = n()) %>%
    rename(Value=!!groups[[i]]) %>%
    mutate(Attr=groups[i])
  Level_Randomization[[i]] <- N
}

Level_Randomization_N <- do.call(rbind,Level_Randomization)
Level_Randomization_Prop <- Level_Randomization_N %>% mutate(Prop=N/762)

Level_Randomization_Label <- data.frame(trimws(gsub("[[:punct:]]","",Level_Randomization_N$Attr),"both"))
names(Level_Randomization_Label) <- "Attribute"

Short.Level <- data.frame(c("Below average","Average","Above average",
                            "Fixed","20%","10%","5%",
                            "Supervisor","Attendance","Changes","Satisfaction",
                            "Rare","Frequent","Moderate",
                            "Low","Average","High",                                
                            "White","Black","Hispanic","Multiracial",
                            "Never","Frequently","Occasionally",
                            "Analysis","Coordination","Interaction","Teamwork"))
names(Short.Level) <- "Short.level"
Level_Randomization_Final <- cbind(Level_Randomization_Prop,Short.Level,Level_Randomization_Label)

ggplot(Level_Randomization_Final, aes(x = Short.level, y = Prop)) +
  geom_col(size=1) +
  facet_wrap(~ Attribute, scales = "free") +
  ylab("Percentage") +
  xlab("Value") +
  theme(axis.text.x = element_text(angle = 8, vjust = 0.5, hjust=0.5))


# 5. Risk Attitude (Appendix Table A8-A9)
# 5.1 Risk Attitude Difference by Race
PFP.Choice.Race.Risk.Test <- P4P_Clean %>% distinct(ID,Race_B,Race_C,GeneralRiskNum,LotteryRisk)
#shapiro.test(PFP.Choice.Race.Risk.Test$GeneralRiskNum)
#shapiro.test(PFP.Choice.Race.Risk.Test$LotteryRisk)
Test1 <- group_by(PFP.Choice.Race.Risk.Test, Race_B) %>%
  summarise(count = n(),
            mean = mean(GeneralRiskNum, na.rm = TRUE),
            sd = sd(GeneralRiskNum, na.rm = TRUE),
            median = median(GeneralRiskNum, na.rm = TRUE),
            IQR = IQR(GeneralRiskNum, na.rm = TRUE))
Test2 <- group_by(PFP.Choice.Race.Risk.Test, Race_B) %>%
  summarise(count = n(),
            mean = mean(LotteryRisk, na.rm = TRUE),
            sd = sd(LotteryRisk, na.rm = TRUE),
            median = median(LotteryRisk, na.rm = TRUE),
            IQR = IQR(LotteryRisk, na.rm = TRUE))
PFP.Choice.Race.Risk.Test.White <- PFP.Choice.Race.Risk.Test %>% subset(Race_B=="White")
PFP.Choice.Race.Risk.Test.Minority <- PFP.Choice.Race.Risk.Test %>% subset(Race_B=="Minority")
wilcox.test(PFP.Choice.Race.Risk.Test.White$GeneralRiskNum,PFP.Choice.Race.Risk.Test.Minority$GeneralRiskNum,paired=FALSE)
wilcox.test(PFP.Choice.Race.Risk.Test.White$LotteryRisk,PFP.Choice.Race.Risk.Test.Minority$LotteryRisk,paired=FALSE)

Test3 <- group_by(PFP.Choice.Race.Risk.Test, Race_C) %>%
  summarise(count = n(),
            mean = mean(GeneralRiskNum, na.rm = TRUE),
            sd = sd(GeneralRiskNum, na.rm = TRUE),
            median = median(GeneralRiskNum, na.rm = TRUE),
            IQR = IQR(GeneralRiskNum, na.rm = TRUE))
Test4 <- group_by(PFP.Choice.Race.Risk.Test, Race_C) %>%
  summarise(count = n(),
            mean = mean(LotteryRisk, na.rm = TRUE),
            sd = sd(LotteryRisk, na.rm = TRUE),
            median = median(LotteryRisk, na.rm = TRUE),
            IQR = IQR(LotteryRisk, na.rm = TRUE))
PFP.Choice.Race.Risk.Test.noAnswer <- PFP.Choice.Race.Risk.Test %>% subset(Race_B!="Prefer not to answer")
kruskal.test(GeneralRiskNum ~ Race_C, data = PFP.Choice.Race.Risk.Test.noAnswer)
kruskal.test(LotteryRisk ~ Race_C, data = PFP.Choice.Race.Risk.Test.noAnswer)


# 5.2 Risk Attitude Difference by Gender
PFP.Choice.Gender.Risk.Test <- P4P_Clean %>% subset(Gender=="Female" | Gender=="Male") %>% 
  distinct(ID,Gender,GeneralRiskNum,LotteryRisk)
Test5 <- group_by(PFP.Choice.Gender.Risk.Test, Gender) %>%
  summarise(count = n(),
            mean = mean(GeneralRiskNum, na.rm = TRUE),
            sd = sd(GeneralRiskNum, na.rm = TRUE),
            median = median(GeneralRiskNum, na.rm = TRUE),
            IQR = IQR(GeneralRiskNum, na.rm = TRUE))
Test6 <- group_by(PFP.Choice.Gender.Risk.Test, Gender) %>%
  summarise(count = n(),
            mean = mean(LotteryRisk, na.rm = TRUE),
            sd = sd(LotteryRisk, na.rm = TRUE),
            median = median(LotteryRisk, na.rm = TRUE),
            IQR = IQR(LotteryRisk, na.rm = TRUE))
PFP.Choice.Gender.Risk.Test.Male <- PFP.Choice.Gender.Risk.Test %>% subset(Gender=="Male")
PFP.Choice.Gender.Risk.Test.Female <- PFP.Choice.Gender.Risk.Test %>% subset(Gender=="Female")
wilcox.test(PFP.Choice.Gender.Risk.Test.Male$GeneralRiskNum,PFP.Choice.Gender.Risk.Test.Female$GeneralRiskNum,paired=FALSE)
wilcox.test(PFP.Choice.Gender.Risk.Test.Male$LotteryRisk,PFP.Choice.Gender.Risk.Test.Female$LotteryRisk,paired=FALSE)


# 5.3 Risk Attitude Difference by Age
PFP.Choice.Age.Risk.Test <- P4P_Clean %>% distinct(ID,Age,Age_B,GeneralRiskNum,LotteryRisk)
Test7 <- group_by(PFP.Choice.Age.Risk.Test, Age_B) %>%
  summarise(count = n(),
            mean = mean(GeneralRiskNum, na.rm = TRUE),
            sd = sd(GeneralRiskNum, na.rm = TRUE),
            median = median(GeneralRiskNum, na.rm = TRUE),
            IQR = IQR(GeneralRiskNum, na.rm = TRUE))
Test8 <- group_by(PFP.Choice.Age.Risk.Test, Age_B) %>%
  summarise(count = n(),
            mean = mean(LotteryRisk, na.rm = TRUE),
            sd = sd(LotteryRisk, na.rm = TRUE),
            median = median(LotteryRisk, na.rm = TRUE),
            IQR = IQR(LotteryRisk, na.rm = TRUE))
PFP.Choice.Age.Risk.Test.Young <- PFP.Choice.Age.Risk.Test %>% subset(Age_B==0)
PFP.Choice.Age.Risk.Test.Old <- PFP.Choice.Age.Risk.Test %>% subset(Age_B==1)
wilcox.test(PFP.Choice.Age.Risk.Test.Young$GeneralRiskNum,PFP.Choice.Age.Risk.Test.Old$GeneralRiskNum,paired=FALSE)
wilcox.test(PFP.Choice.Age.Risk.Test.Young$LotteryRisk,PFP.Choice.Age.Risk.Test.Old$LotteryRisk,paired=FALSE)


# 6. Potential Job Seekers (Student & Unemployed) (Appendix Table A10-A14)
table(P4P_Clean$Employment,exclude=NULL)
summary(P4P_Clean$Age)

P4P_Data_Seeker <- P4P_Clean %>%
  subset(Employment=="Student" | Employment=="Unemployed" | Employment=="Employed") %>%
  subset(Age<=50)
summary(P4P_Data_Seeker$Age)

JS1 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
            `Current community involvement` + `Community income` + `Community demographics` +
            `Overtime work` + `Key job task`, 
          data=P4P_Data_Seeker)
summary(JS1)
JS1.cl <- cluster.vcov(JS1, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JS1, JS1.cl)
JS1.Twoway.se <- sqrt(diag(JS1.cl))


# 6.1 Levels of P4P
JS2 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
            `Current community involvement` + `Community income` + `Community demographics` +
            `Overtime work` + `Key job task`, 
          data=P4P_Data_Seeker)
summary(JS2)
JS2.cl <- cluster.vcov(JS2, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JS2, JS2.cl)
JS2.Twoway.se.NotCollapse <- sqrt(diag(JS2.cl))


# 6.2 Interaction w/ Goal based evaluation
JS3 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
            `Current community involvement` + `Community income` + `Community demographics` +
            `Overtime work` + `Key job task` + `Performance bonuses binary`:`Goal based evaluation`, 
          data=P4P_Data_Seeker)
summary(JS3)
JS3.cl <- cluster.vcov(JS3, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JS3, JS3.cl)
JS3.Twoway.se.PG <- sqrt(diag(JS3.cl))


# 6.3 Interaction w/ Risk Attitude
JS4 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
            `Current community involvement` + `Community income` + `Community demographics` +
            `Overtime work` + `Key job task` + GeneralRiskB + `Performance bonuses binary`:GeneralRiskB, 
          data=P4P_Data_Seeker)
summary(JS4)
JS4.cl <- cluster.vcov(JS4, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JS4, JS4.cl)
JS4.Twoway.se.PRG <- sqrt(diag(JS4.cl))

JS5 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
            `Current community involvement` + `Community income` + `Community demographics` +
            `Overtime work` + `Key job task` + LotteryRiskB + `Performance bonuses binary`:LotteryRiskB, 
          data=P4P_Data_Seeker)
summary(JS5)
JS5.cl <- cluster.vcov(JS5, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JS5, JS5.cl)
JS5.Twoway.se.PRL <- sqrt(diag(JS5.cl))

## Report (Table A10)
stargazer(JS1,JS2,JS3,JS4,JS5,
          se=list(JS1.Twoway.se,JS2.Twoway.se.NotCollapse,JS3.Twoway.se.PG,JS4.Twoway.se.PRG,JS5.Twoway.se.PRL), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%)",
                             "Performance bonuses: moderate (10%)",
                             "Performance bonuses: small (5%)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "(A) X (B)",
                             "General Risk: High Risk", 
                             "(A) X General Risk: High Risk",
                             "Lottery Risk: High Risk",
                             "(A) X Lottery Risk: High Risk"),
          out="Seeker_P4P.txt") 


# 6.4 Race 
table(P4P_Data_Seeker$Race_B,exclude=NULL)
table(P4P_Data_Seeker$Race_C,exclude=NULL)

P4P_Data_Seeker_Race <- P4P_Data_Seeker %>%
  subset(Race_B!="Prefer not to answer")

JSR1 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Race_B + `Performance bonuses binary`*Race_B, 
           data=P4P_Data_Seeker_Race)
summary(JSR1)
JSR1.cl <- cluster.vcov(JSR1, cbind(P4P_Data_Seeker_Race$ID,P4P_Data_Seeker_Race$Treatment.ID))
coeftest(JSR1, JSR1.cl)
JSR1.Twoway.se.RaceB <- sqrt(diag(JSR1.cl))


# 6.4.1 Levels of P4P
JSR2 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Race_B + `Performance bonuses`:Race_B, 
           data=P4P_Data_Seeker_Race)
summary(JSR2)
JSR2.cl <- cluster.vcov(JSR2, cbind(P4P_Data_Seeker_Race$ID, P4P_Data_Seeker_Race$Treatment.ID))
coeftest(JSR2, JSR2.cl)
JSR2.Twoway.se.RaceB.Prop <- sqrt(diag(JSR2.cl))

# 6.4.2 Interaction w/ Goal based evaluation
JSR3 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Race_B + 
             `Performance bonuses binary`:Race_B + `Goal based evaluation`:Race_B +
             `Performance bonuses binary`:`Goal based evaluation` +
             `Performance bonuses binary`:`Goal based evaluation`:Race_B, 
           data=P4P_Data_Seeker_Race)
summary(JSR3)
JSR3.cl <- cluster.vcov(JSR3, cbind(P4P_Data_Seeker_Race$ID, P4P_Data_Seeker_Race$Treatment.ID))
coeftest(JSR3, JSR3.cl)
JSR3.Twoway.se.EVAL <- sqrt(diag(JSR3.cl))

# 6.4.3 Interaction w/ Risk Attitude
JSR4 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Race_B + GeneralRiskB +
             `Performance bonuses binary`*Race_B + `Performance bonuses binary`*GeneralRiskB, 
           data=P4P_Data_Seeker_Race)
summary(JSR4)
JSR4.cl <- cluster.vcov(JSR4, cbind(P4P_Data_Seeker_Race$ID, P4P_Data_Seeker_Race$Treatment.ID))
coeftest(JSR4, JSR4.cl)
JSR4.Twoway.se.RaceB.RiskG <- sqrt(diag(JSR4.cl))

JSR5 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Race_B + LotteryRiskB +
             `Performance bonuses binary`*Race_B + `Performance bonuses binary`*LotteryRiskB, 
           data=P4P_Data_Seeker_Race)
summary(JSR5)
JSR5.cl <- cluster.vcov(JSR5, cbind(P4P_Data_Seeker_Race$ID, P4P_Data_Seeker_Race$Treatment.ID))
coeftest(JSR5, JSR5.cl)
JSR5.Twoway.se.RaceB.RiskL <- sqrt(diag(JSR5.cl))


## Report (Table A11)
stargazer(JSR1,JSR2,JSR3,JSR4,JSR5,
          se=list(JSR1.Twoway.se.RaceB,JSR2.Twoway.se.RaceB.Prop,JSR3.Twoway.se.EVAL,JSR4.Twoway.se.RaceB.RiskG,JSR5.Twoway.se.RaceB.RiskL), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%) (A1)",
                             "Performance bonuses: moderate (10%) (A2)",
                             "Performance bonuses: small (5%) (A3)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "Race - binary: Minority (C)",
                             "General Risk: High Risk",
                             "Lottery Risk: High Risk",
                             "(A) X (C)",
                             "(A1) X (C)",
                             "(A2) X (C)",
                             "(A3) X (C)",
                             "(B) X (C)",
                             "(A) X (B)",
                             "(A) X (B) X (C)",
                             "(A) X General Risk: High Risk",
                             "(A) X Lottery Risk: High Risk"),
          out="Seeker_Race_P4P.txt") 


# 6.5 Race - Specific
JSR6 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Race_C + `Performance bonuses binary`:Race_C, 
           data=P4P_Data_Seeker_Race)
summary(JSR6)
JSR6.cl <- cluster.vcov(JSR6, cbind(P4P_Data_Seeker_Race$ID, P4P_Data_Seeker_Race$Treatment.ID))
coeftest(JSR6, JSR6.cl)
JSR6.Twoway.se.RaceC <- sqrt(diag(JSR6.cl))

JSR7 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Race_C + `Performance bonuses`:Race_C, 
           data=P4P_Data_Seeker_Race)
summary(JSR7)
JSR7.cl <- cluster.vcov(JSR7, cbind(P4P_Data_Seeker_Race$ID, P4P_Data_Seeker_Race$Treatment.ID))
coeftest(JSR7, JSR7.cl)
JSR7.Twoway.se.RaceC.Prop <- sqrt(diag(JSR7.cl))

## Report (Table A12)
stargazer(JSR6, JSR7,
          se=list(JSR6.Twoway.se.RaceC,JSR7.Twoway.se.RaceC.Prop), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%) (A1)",
                             "Performance bonuses: moderate (10%) (A2)",
                             "Performance bonuses: small (5%) (A3)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "Race - Asian", "Race - Black", "Race - Hispanic", "Race - Other",
                             "(A) X Race - Asian", "(A) X Race - Black", "(A) X Race - Hispanic", "(A) X Race - Other",
                             "(A1) X Race - Asian", "(A2) X Race - Asian", "(A3) X Race - Asian",
                             "(A1) X Race - Black", "(A2) X Race - Black", "(A3) X Race - Black",
                             "(A1) X Race - Hispanic", "(A2) X Race - Hispanic", "(A3) X Race - Hispanic",
                             "(A1) X Race - Other", "(A2) X Race - Other", "(A3) X Race - Other"),
          out="Seeker_Race_Separate_P4P.txt") 


# 6.6. Gender
P4P_Data_Seeker_Gender <- P4P_Data_Seeker %>%
  subset(Gender=="Female" | Gender=="Male") %>%
  mutate(Gender=factor(Gender))

JSG1 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + `Performance bonuses binary`:Gender, 
           data=P4P_Data_Seeker_Gender)
summary(JSG1)
JSG1.cl <- cluster.vcov(JSG1, cbind(P4P_Data_Seeker_Gender$ID, P4P_Data_Seeker_Gender$Treatment.ID))
coeftest(JSG1, JSG1.cl)
JSG1.Twoway.se.Gender <- sqrt(diag(JSG1.cl))

# 6.6.1 Levels of P4P
JSG2 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + `Performance bonuses`:Gender, 
           data=P4P_Data_Seeker_Gender)
summary(JSG2)
JSG2.cl <- cluster.vcov(JSG2, cbind(P4P_Data_Seeker_Gender$ID, P4P_Data_Seeker_Gender$Treatment.ID))
coeftest(JSG2, JSG2.cl)
JSG2.Twoway.se.Gender.Prop <- sqrt(diag(JSG2.cl))

# 6.6.2 Interaction w/ Goal based evaluation
JSG3 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + `Performance bonuses binary`*`Goal based evaluation`*Gender, 
           data=P4P_Data_Seeker_Gender)
summary(JSG3)
JSG3.cl <- cluster.vcov(JSG3, cbind(P4P_Data_Seeker_Gender$ID, P4P_Data_Seeker_Gender$Treatment.ID))
coeftest(JSG3, JSG3.cl)
JSG3.Twoway.se.EVAL.Gender <- sqrt(diag(JSG3.cl))

# 6.6.3 Interaction w/ Risk Attitude
JSG4 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + GeneralRiskB +
             `Performance bonuses binary`:Gender + `Performance bonuses binary`:GeneralRiskB, 
           data=P4P_Data_Seeker_Gender)
summary(JSG4)
JSG4.cl <- cluster.vcov(JSG4, cbind(P4P_Data_Seeker_Gender$ID, P4P_Data_Seeker_Gender$Treatment.ID))
coeftest(JSG4, JSG4.cl)
JSG4.Twoway.se.Gender.RiskG <- sqrt(diag(JSG4.cl))

JSG5 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Gender + LotteryRiskB +
             `Performance bonuses binary`:Gender + `Performance bonuses binary`:LotteryRiskB, 
           data=P4P_Data_Seeker_Gender)
summary(JSG5)
JSG5.cl <- cluster.vcov(JSG5, cbind(P4P_Data_Seeker_Gender$ID, P4P_Data_Seeker_Gender$Treatment.ID))
coeftest(JSG5, JSG5.cl)
JSG5.Twoway.se.Gender.RiskL <- sqrt(diag(JSG5.cl))


## Report (Table A13)
stargazer(JSG1,JSG2,JSG3,JSG4,JSG5,
          se=list(JSG1.Twoway.se.Gender,JSG2.Twoway.se.Gender.Prop,JSG3.Twoway.se.EVAL.Gender,JSG4.Twoway.se.Gender.RiskG,JSG5.Twoway.se.Gender.RiskL), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%) (A1)",
                             "Performance bonuses: moderate (10%) (A2)",
                             "Performance bonuses: small (5%) (A3)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "Gender: Women (C)",
                             "(A) X (B)",
                             "General Risk: High Risk", "Lottery Risk: High Risk",
                             "(A) X (C)",
                             "(A1) X (C)",
                             "(A2) X (C)",
                             "(A3) X (C)",
                             "(B) X (C)",
                             "(A) X (B) X (C)",
                             "(A) X General Risk: High Risk", "(A) X Lottery Risk: High Risk"),
          out="Seeker_Gender_P4P.txt") 


# 6.7 Age
JSA1 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + `Performance bonuses binary`:Age_B, 
           data=P4P_Data_Seeker)
summary(JSA1)
JSA1.cl <- cluster.vcov(JSA1, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JSA1, JSA1.cl)
JSA1.Twoway.se.Age <- sqrt(diag(JSA1.cl))

# 6.7.1 Levels of P4P
JSA2 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + `Performance bonuses`:Age_B, 
           data=P4P_Data_Seeker)
summary(JSA2)
JSA2.cl <- cluster.vcov(JSA2, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JSA2, JSA2.cl)
JSA2.Twoway.se.Age.Prop <- sqrt(diag(JSA2.cl))

# 6.7.2 Interaction w/ Goal based evaluation
JSA3 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + `Performance bonuses binary`*`Goal based evaluation`*Age_B, 
           data=P4P_Data_Seeker)
summary(JSA3)
JSA3.cl <- cluster.vcov(JSA3, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JSA3, JSA3.cl)
JSA3.Twoway.se.EVAL.Age <- sqrt(diag(JSA3.cl))

# 6.7.3 Interaction w/ Risk Attitude
JSA4 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + GeneralRiskB +
             `Performance bonuses binary`:Age_B + `Performance bonuses binary`:GeneralRiskB, 
           data=P4P_Data_Seeker)
summary(JSA4)
JSA4.cl <- cluster.vcov(JSA4, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JSA4, JSA4.cl)
JSA4.Twoway.se.Age.RiskG <- sqrt(diag(JSA4.cl))

JSA5 <- lm(Chosen_Job ~ `Total pay` + `Performance bonuses binary` + `Goal based evaluation` +
             `Current community involvement` + `Community income` + `Community demographics` +
             `Overtime work` + `Key job task` + Age_B + LotteryRiskB +
             `Performance bonuses binary`:Age_B + `Performance bonuses binary`:LotteryRiskB, 
           data=P4P_Data_Seeker)
summary(JSA5)
JSA5.cl <- cluster.vcov(JSA5, cbind(P4P_Data_Seeker$ID, P4P_Data_Seeker$Treatment.ID))
coeftest(JSA5, JSA5.cl)
JSA5.Twoway.se.Age.RiskL <- sqrt(diag(JSA5.cl))

## Report (Table A14)
stargazer(JSA1,JSA2,JSA3,JSA4,JSA5,
          se=list(JSA1.Twoway.se.Age,JSA2.Twoway.se.Age.Prop,JSA3.Twoway.se.EVAL.Age,JSA4.Twoway.se.Age.RiskG,JSA5.Twoway.se.Age.RiskL), type="text", 
          keep.stat=c("n","adj.rsq"), no.space=TRUE, 
          dep.var.caption="", dep.var.labels.include=FALSE,
          covariate.labels=c("Total pay: About average","Total pay: Slightly above average",
                             "Performance bonuses binary: Yes (A)",
                             "Performance bonuses: large (20%) (A1)",
                             "Performance bonuses: moderate (10%) (A2)",
                             "Performance bonuses: small (5%) (A3)",
                             "Goal based evaluation: Yes (B)",
                             "Current community involvement: Frequent participation",
                             "Current community involvement: Moderate participation",
                             "Community income: Average income", "Community income: High income",
                             "Community demographics: Mostly African American",
                             "Community demographics: Mostly Hispanic",
                             "Community demographics: Multiracial",
                             "Overtime work: Frequently required","Overtime work: Occasionally required",
                             "Key job task: Coordination with community groups",
                             "Key job task: Direct interaction with community residents",
                             "Key job task: Teamwork with peers and supervisors",
                             "Age: Older (C)",
                             "(A) X (B)",
                             "General Risk: High Risk", "Lottery Risk: High Risk",
                             "(A) X (C)",
                             "(A1) X (C)",
                             "(A2) X (C)",
                             "(A3) X (C)",
                             "(B) X (C)",
                             "(A) X (B) X (C)",
                             "(A) X General Risk: High Risk", "(A) X Lottery Risk: High Risk"),
          out="Seeker_Age_P4P.txt") 

